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Abstract 

Violent relaxation process of spherical stellar systems is examined by numer- 
ical simulations of shell model. The collapse of uniform density sphere both 
with and without external force is investigated. It is found that time vari- 
ation of mean gravitational potential field induced by external force makes 
the system closer to Lynden-Bell distribution. Our result imply that sys- 
tems which are gravitationary affected by other systems may approach the 
Lynden-Bell distribution. 

Key words: Collisionless stellar dynamics; Elliptical galaxies; Violent 
relaxation 



1 Introduction 

Violent relaxation is a mechanism proposed to explain a similarity and smooth- 
ness of the observed light distribution of elliptical galaxies. This relaxation 
process is caused by explicitly time-dependant mean gravitational potential 
field. Its typical timescale is much shorter than two-body relaxation time 
and comparable with dynamical time. 

According to Lynden-Bell (1967), it is expected that the coarse-grained 
entropy of systems becomes maximum through violent relaxation in ideal 
cases and that Lynden-Bell distribution can be established. The Lynden- 
Bell distribution is given by 

if^ exp{-/?(g-/i)} 

^^'^=^l + exp{-/5(e-^)}' ^'^ 

where f{e) is coarse-grained distribution function, e is energy per unit mass, 
and 1], j3, n are constants. In almost all cases, degeneracy can be neglected 
(/(e) ^ 1]), so Lynden-Bell distribution is approximated by the Maxwell- 
Boltzmann distribution: 

/» = Aexp(-/3£), (2) 

where A = ?7exp(/9/i) = constant. However, in numerical simulations so far 
of one dimensional gravitational N-body systems (sheet systems), they did 
not reach Lynden-bell distribution through violent relaxation ( Cuperman 
et al. 1969; Hohl and Campbell 1968; Goldstein et al. 1969; Yamashiro 



et al. 1992). It remains unclear why Lynden-Bell distribution cannot be 
obtained in these simulations. Moreover it has not been fully investigated 
what distributions are obtained in other gravitational N-body systems (shell 
systems, particle systems) through violent relaxation. 

It is widely believed that the Lynden-Bell distribution is not realized in N- 
body simulations (e.g., Funato et al. 1992a, 1992b). One of the reasons is that 
the fluctuations in the gravitational potential decays before the relaxation 
process is completed. To investigate if this is true, we add artificial time 
(and space) varying potentials and see if the Lynden-Bell distribution is 
established. As a first step, we adopt simple spherical symmetric models 
since spherical symmetric models are more realistic than sheet models. More 
realistic models will be studied in subsequent papers. 

2 Numerical Method 

In our simulations, shell models (Henon 1964) are adapted. In spherical sym- 
metric stellar systems, distribution function can be expressed as f{r,u,v,t), 
where r the radial distance, u the radial velocity, v the tangential velocity, 
and t is the time, respectively. So distribution function is represented by A^ 
points in (r, u, v) space and each point evolves according to the equations of 
motion. 

^ = «^' (3) 



dm ^ A\ GMj d(j)eAr,t) 
dt rf rf dr '' 

where G the gravitational constant, Ai = ViVi = constant (with respect to 
time) is the angular momentum of each shell, Mi the mass interior to r,, and 
4'ex{i^,t) is the potential of external force. 

We set G = M(total mass)= 1 and A^ = 10000 . And each shell has 
the same mass (= 1/A^) . The integration is carried out using a leap-frog 
method. 

When we set (pexif^jt) = 0, we set integration time-step small enough for 
the maximum error of total energy to be of the order of 0.01% after lOOtcr 
. tcr is the crossing time which is typical timescale for shells to cross typical 
radius of the system. In our simulations, tcr ~ 1 in this system of units. 
When (pexij'it) 7^ 0, we set integration time-step same as the model whose 
initial condition is common. 

3 Models 

We calculate collapses of uniform density sphere whose initial virial ratio is 
0.5 and 0.1. In some models external gravitational potential field are added 
after systems are settled down to study equilibrium states. 

Initial densities and velocity distributions are common in all models. 

• p{r) = const. (0 < r < 2) . 

• Both u and v obey isotropic Maxwell distribution. 



For simplicity, we assume that both time and spatial dependence of ex- 
ternal potential field is sinusoidal. We set (pex as follows, 

rAsin(27rt/P)sin(27rr/A) (20 < t < 20 + 5P) 

In (^ P and A are parameters. We set external potential to survive 
for much shorter time than two-body relaxation time which is comparable 
to Ntcr ( see appendix 1 ). We set typical intensity of external potential 
comparable to that of inertial potential. In our simulations, we set A=0.5 
when initial virial ratio is 0.5, and A=1.36 when initial virial ratio 0.1. 

In table 1 the parameters for each models are listed. In this table both P 
and A are normalized in the units of R and i, which are typical lengthscale 
and time scale of the post collapsed system. In models whose initial virial 
ratio is 0.5, R = 1.34 and i = 1.53, and in models whose initial virial ratio is 
0.1, R = 0.494 and i = 0.345 ( see appendix 2 ). Models are summarized in 
table 1. 

4 Results 

In models A-0 and B-0, the system approached stationary states within a 
few crossing times and variation of virial ratio and variation of distribution 
function was damped(figure 1 and figure 2). In figure 2, the mean value of the 
virial ratio at the stationary stages is greater than 1.0 (about 1.35) because 
the existence of shells which have positive energies cannot be negligible. 
Figure 1 Figure 2 



In other models, the systems also approached stationary states within a 
few crossing time after variation of external potential stopped. In table 2 
ratio of mass whose energy is more than zero in stationary states is listed in 
each model. 

In figures 3, figures 4 and figures 5, distribution function of stationary 
states of model A-0, model B-0 and model A-1 are shown by the solid lines. 
We derive distribution function from the simulation data as follows. We 
divide energy space between e = Smin and e = into 100 bins, where Smin 
is the minimum value of the potential energy of the system. Therefor the 
i-th bin corresponds to the energy range between Si — Ae/2 and Si + A£:/2, 
where Si = Emm + (^ ~ 1/2)A£ and Ae = — £mm/100. Then we count the 
number of shells contained the i-th bin and it is divided by the phase space 
volume which corresponds to the energy range to get the /(£«). Anisotropy 
of velocity distribution cannot harm this procedure because energy of each 
shell doesn't depend on the direction of velocity but only on the absolute 
value. 



Figure 3 Figure 4 Figure 5 



We estimate the deviation from Lynden-Bell distribution of each model as 
follows. Unfortunately, spherical stellar systems with Maxwell distribution 
(isothermal distribution) have infinite mass. So the distribution expected 
after violent relaxation in real spherical stellar systems is not Maxwellian 
itself, but that of the King model which has 'cut-off' at higher energy end. 



So distribution functions in stationary states are fitted by the form, 

/» = pi(27ra^)-^/^{exp(^i±^) - 1}, (6) 

where, pi, a, $0 are fitting parameters. In figure 3, figure 4 and figure 5, 
distribution functions of the best fit King-models are showed by the short 
dash hues. 

Then we define D^ as follows and use it for the estimation of the deviation 
from Lynden-Bell distribution. 



Nh 



^^- E-0 E 



-1.^ fdatai^i) — ffiti^i 



n. 

i=l i=l 



(7) 



ffitiSi) 

where, fdata is the distribution function from simulation data, ffu is the best 
fit function of the data in the form of (|^), Nf, is the number of bins (in this 
paper A^^ = 100) and n^ is the number of shells contained in the i-th bin. 
Fluctuations from the finiteness of shell number are proportional to square 
root of shell number. Therefore, the summation is carried out weighted with 
rii in the numerator of equation (0) . 

In table 3, fitting parameters ( pi, a^, ^q ), c, which is the concentration 
parameter of the corresponding King model obtained from the fitting result 
and D^ of each model are listed. Compared between Model A-0 and B-0, 
both of which are without external force, D^ of B-0 is smaller than that of 
A-0. Initial condition of model B-0 is farther from dynamical equilibrium 
state than that of A-0, so time variation of gravitational potential is larger 
and relaxation is more effective. In all models whose initial virial ratios are 



0.5 with external force (from A-1 to A-9), D^ becomes smaller than that of 
model A-0 which has the same initial condition but without external force. 
The similar tendency can be seen in almost models whose initial virial ratio 
is 0.1 (from B-0 to B-1) except B-4 and B-7. 

5 Discussion 

We examined the violent relaxation process of spherical stellar systems and 
the influence on it of time variation of external gravitational potential fields 
by using shell models. It is found that time variation of the mean potential 
field mainly caused by external gravitational field makes distribution function 
closer to Maxwellian. Especially in some models (Al, A4, A7, B8, B9) 
the deviations from the King model become quite small. Though it is not 
clear why the deviations are small in these models, it shows that the violent 
relaxation is quite efficient under some conditions. 

So, if time variation of mean gravitational potential fields continues for a 
long enough time, systems can reach Lynden-Bell distribution. 

It is rare that no gravitational infiuence of other galaxies is exercised on 
galaxies. Especially in galaxies in clusters, this effect plays an important role 
on their dynamical evolution. Though it is not clear the real gravitational 
perturbations are how much effective to cause the violent relaxation, there 
are possibilities that the distribution function of the real elliptical galaxies 
become close to the King models by the external forces. 



9 



We would like to thank T.Hayashi for his contribution to developing nu- 
merical code and M. Shimada for helpful discussions. Numerical computa- 
tions in this work were carried out with workstations at Yukawa Institute 
for Theoretical Physics. This work is in part supported by Research Fellow- 
ships of the Japan Society for the Promotion of Science for Young Scientists 
(M.T.). 

Appendix 1. Estimation of Two-Body Re- 
laxation Time of Shell Systems 

We assume that there are N shells which have the same mass of m . Then 
the mean energy of each shell e is as follows from the Virial theorem. 

(Al.l) 



R ' 

where G is the gravitational constant and ^ is a typical radius of the system. 
And the mean energy change of each shell in one shell-crossing-event is 

\Ae\ ~ -^. (A1.2) 

Relaxation process can be regarded as one-dimensional random walks in 
the energy space whose one step-size is equal to |Ae| . Since the system can 
be regarded as relaxed when the variance of displacement in the energy space 
becomes comparable to the square of typical energy of each shell, number of 
crossing-event required for relaxation, n^e/, is calculated as follows. 

e^r^nrei\M\ (A1.3) 
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Urel ~ {^Y ~ N\ (A1.4) 

Then relaxation time of shell systems, tj-eu is 

f R 

trel ~ nrel- ~ A^— ~ Ntcr, (Al.5) 

W V 

where f = ^/A^ is the mean separation between shells and v is the mean 
velocity of shells. 

Appendix 2. Changes of Length and Time Scale 

through Collapse 

We assume that the initial total mass is M, the initial typical length scale is 
-Ro; and that the initial velocity dispersion is Vq. Then total kinetic energy, 
To, total potential energy, Wq, and virial ratio, Vq, of initial states are. 

To = ^Mvl (A2.1) 

GM^ , , , 

W^o = -^, (A2.2) 

ito 



Vi 







2Tf 



GM' 



(A2.3) 



where G is the gravitational constant. 

When the system reached virial equilibrium after the collapse, we let the 
typical length scale i?, velocity dispersion f^, and the total mass M(l — x), 
where x is the escape rate of mass. Then total kinetic energy, T, total 
potential energy, W ^ and virial ratio, V, of the bounded portion of the system 
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are, 

T = -M{l-x)v^, (A2.4) 

If we let the energy which escaped mass possess Ai5, according to energy 
conservation law, 



Tq + Wq = T + W + AE, (A2.7) 



then, 

v^ = {2-Vro){l-x)-\l + -—) (A2.8) 

o A F 

Vro \Eo\ 

AE 
R = Roi2-Vro)-\l-xy{l + --)-\ (A2.10) 

t = ^L(— ^)2/3(i_a;f/2(i + ^)-3/2 (A2.11) 

- . Ylil n_a;)5/2(i + ^)-3/2^ (A2.12) 



where |ii^o| = G'M^(1 — \/ro/2)/i?o is the absolute value of initial total energy 
of the system, t = R/v is the typical timescale of the system after the system 
has collapsed, and to = Rq/vq is that of initial state. 

In our simulation, when Vq = 0.5, x and AE/\Eq\ is much smaller than 
unity and can be negligible. On the other hand, when Vq = 0.1, a; = 0.32 and 
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/\E/\Eq\ = 0.45 . With our unit { G = M = 1 ) and the initial conditions 
described in section 3, when Vq = 0.5, R = 1.34 and t = 1.53 from ( |A2.10| ) 
and ( [A2Tl| ). When Vq = 0.1, R = 0.494 and t = 0.345. 
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Table 1. Parameters of Models 



Model 


Initial V.R. 


External P.F. 


P/i 


\/R 


A-0 


0.5 


No 


- 


- 


A-1 


0.5 


Yes 


0.5 


0.5 


A-2 


0.5 


Yes 


0.5 


1.0 


A-3 


0.5 


Yes 


0.5 


2.0 


A-4 


0.5 


Yes 


1.0 


0.5 


A-5 


0.5 


Yes 


1.0 


1.0 


A-6 


0.5 


Yes 


1.0 


2.0 


A-7 


0.5 


Yes 


2.0 


0.5 


A-8 


0.5 


Yes 


2.0 


1.0 


A-9 


0.5 


Yes 


2.0 


2.0 


B-0 


0.1 


No 


- 


- 


B-1 


0.1 


Yes 


0.5 


0.5 


B-2 


0.1 


Yes 


0.5 


1.0 


B-3 


0.1 


Yes 


0.5 


2.0 


B-4 


0.1 


Yes 


1.0 


0.5 


B-5 


0.1 


Yes 


1.0 


1.0 


B-6 


0.1 


Yes 


1.0 


2.0 


B-7 


0.1 


Yes 


2.0 


0.5 


B-8 


0.1 


Yes 


2.0 


1.0 


B-9 


0.1 


Yes 


2.0 


2.0 



15 



Table 2. Ratio of mass whose energy is more than zero in stationary 
states. 



Model 


Rate 


Model 


Rate 


A-0 


0.0003 


B-0 


0.33 


A-1 


0.78 


B-1 


0.85 


A-2 


0.11 


B-2 


0.39 


A-3 


0.0029 


B-3 


0.31 


A-4 


0.41 


B-4 


0.60 


A-5 


0.78 


B-5 


0.82 


A-6 


0.20 


B-6 


0.43 


A-7 


0.18 


B-7 


0.42 


A-8 


0.53 


B-8 


0.64 


A-9 


0.75 


B-9 


0.80 
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Table 3. Fitting parameters ( pi, a^, $0 ) and concentration parameter, 
c, and D^ of each model 



Model 


Pi 




a' 




$0 




c 


D' 


A-0 


2.43 X 


10-4 


0.181 




-2.25 


X 


10-2 


1.52 


1.26 


A-1 


6.31 X 


10-9 


9.80 X 10- 


-3 


-1.56 


X 


10-3 


2.44 


0.360 


A-2 


6.05 X 


10-5 


0.102 




-1.29 


X 


10-2 


1.54 


0.725 


A-3 


2.13 X 


10-4 


0.168 




-2.10 


X 


10-2 


1.51 


0.954 


A-4 


4.32 X 


10-^ 


2.71 X 10" 


-2 


-4.18 


X 


10-3 


2.05 


0.427 


A-5 


5.71 X 


10-9 


9.80 X 10- 


-3 


-1.67 


X 


10-3 


2.33 


0.588 


A-6 


2.57 X 


10-6 


4.27 X 10- 


-2 


-6.15 


X 


10-3 


1.86 


1.19 


A-7 


1.60 X 


10-5 


7.40 X 10- 


-2 


-8.95 


X 


10-3 


1.67 


0.363 


A-8 


4.81 X 


10-6 


5.32 X 10- 


-2 


-6.62 


X 


10-3 


1.61 


0.506 


A-9 


2.22 X 


10-s 


1.33 X 10- 


-2 


-2.11 


X 


10-3 


2.13 


0.852 


B-0 


1.80 X 


10-2 


0.541 




-6.48 


X 


10-2 


1.46 


0.923 


B-1 


1.20 X 


10-9 


7.19 X 10" 


-3 


-1.52 


X 


10-3 


3.01 


0.633 


B-2 


2.12 X 


10-3 


0.227 




-2.67 


X 


10-2 


1.45 


0.720 


B-3 


1.50 X 


10-2 


0.485 




-5.60 


X 


10-2 


1.41 


0.726 


B-4 


1.55 X 


10-5 


5.32 X 10- 


-2 


-8.65 


X 


10-3 


2.11 


1.04 


B-5 


9.33 X 


10-9 


1.04 X 10- 


-2 


-2.31 


X 


10-3 


3.01 


0.534 


B-6 


6.17 X 


10-4 


0.149 




-1.82 


X 


10-2 


1.49 


0.779 


B-7 


4.93 X 


10-4 


0.151 




-1.81 


X 


10-2 


1.67 


1.29 


B-8 


9.48 X 


10-^ 


3.32 X 10- 


-2 


-6.35 


X 


10-3 


2.60 


0.387 


B-9 


2.12 X 


10-s 


1.21 X 10- 


-2 


-2.48 


X 


10-3 


2.78 


0.250 
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Figure captions 

Fig. 1. Time variation of virial ratio of model A-0. Tlie system approaclied 
stationary state witliin a few crossing times and variation of virial ratio 
was damped. 

Fig. 2. Same as figure 1, but for model B-0. the mean of virial ratio of 
stationary state is greater than 1.0 (about 1.35) because the existence 
of shells which have energy of more than zero cannot be negligible. 

Fig. 3. Distribution function of stationary states of model A-0 (by the solid 
line) and that of the best fit King-model (by the short dash line). 

Fig. 4. Same as figure 3, but for model B-0. 

Fig. 5. Same as figure 3, but for model A-1. f{e) obeys Maxwellian. The 
deviation from Maxwellian of lower energy region is due to the statisti- 
cal fiuctuation because number of shells in this energy range is rather 
small. 
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